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Abstract. Network theory provides various tools for investigating the structural or functional topology of 
many complex systems found in nature, technology and society. Nevertheless, it has recently been realised 
that a considerable number of systems of interest should be treated, more appropriately, as interacting 
networks or networks of networks. Here we introduce a novel graph-theoretical framework for studying the 
interaction structure between subnetworks embedded within a complex network of networks. This frame- 
work allows us to quantify the structural role of single vertices or whole subnetworks with respect to the 
interaction of a pair of subnetworks on local, mesoscopic and global topological scales. 
Climate networks have recently been shown to be a powerful tool for the analysis of climatological data. 
Applying the general framework for studying interacting networks, we introduce coupled climate subnet- 
works to represent and investigate the topology of statistical relationships between the fields of distinct 
climatological variables. Using coupled climate subnetworks to investigate the terrestrial atmosphere's 
three-dimensional geopotential height field uncovers known as well as interesting novel features of the at- 
mosphere's vertical stratification and general circulation. Specifically, the new measure "cross-betweenness" 
identifies regions which are particularly important for mediating vertical wind field interactions. The 
promising results obtained by following the coupled climate subnetwork approach present a first step 
towards an improved understanding of the Earth system and its complex interacting components from a 
network perspective. 



1 Introduction 

Complex networks are recognised as a structurally sim- 
ple, yet powerful representation of the manifold systems 
of intricately interacting elements found in nature, tech- 
nology and human society [h{3]. Drawing on ideas from 
mathematical graph theory and statistical physics, com- 
plex network theory allows a detailed and quantitative 
investigation of the interaction topology of networked sys- 
tems Hi as well as exploring the interplay between network 
structure and dynamics on the interacting elements 15]. 

Most studies have so far concentrated on networks 
where vertices represent single elements or subsystems, 
and edges indicate interactions or relationships between 
vertices. However, it has recently been realised that a con- 
siderably large class of systems of interest warrants a more 
natural representation as interacting networks or networks 
of networks for an appropriate description of their in- 
teraction structure (Fig. [T). Notable examples are repre- 
sentations of the mammalian cortex, where cortical areas 
form complex subnetworks that are themselves linked via 
a complex network topology [6j[7], systems of interacting 
populations of heterogeneous oscillators (si [9] , or mutu- 



ally interdependent infrastructure networks 10 11 such 



as the power grid and electricity consuming communica- 
tion networks [12] . More generally, and particularly if the 
system's representation as an interacting network is a less 
obvious choice than for the aforementioned examples, net- 
works with a pronounced community structure may be 
viewed as interacting networks, where subnetworks are 
constituted by communities or clusters as identified by 
some community detection algorithm 113]. A related but 
distinct concept is presented by layerednetworks [14- 16 , 
where essentially different sets of edges (layers) are con- 
sidered connecting the same substrate of vertices, e.g., de- 
scribing a two-layered transportation network with roads 
forming a layer of physical connections between locations 
and a superposed layer of virtual connections induced by 
actual pathways of traffic flow on the physical layer. 

Recently, climate networks representing the statisti- 
cal similarity structure of a spatio-temporally resolved cli- 
matological field were successfully employed for revealing 
novel aspects of climate dynamics in reanalysis data sets 
and in the results of global climate models [17 - 24 . These 
findings include among others insights into the"! effect of the 
El-Niiio Southern Oscillation (ENSO) on the global corre- 
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Fig. 1. Systems of interacting networks or networks of 
networks are a natural representation of many systems 
found in nature, technology and society. The partition 
into interacting subnetworks (indicated by vertex sym- 
bols of different shape) is either naturally induced by 
the considered problem, e.g., consider interdependent in- 
frastructures 



10-12 



or cortical areas of the mammalian 
brain [6j[7], or may"be generated by a community detec- 
tion algorithm 13 . Such a system is characterised by de- 



pendencies within subnetworks (internal edges, continuous 
lines) as well as interactions between different subnetworks 
(cross edges, dashes lines). 



lation structure of several climatological observables 22 
24 as well as the detection of a backbone of significantly 
increased matter and energy flow in the global surface air 
temperature field 17p8 . Now understanding the complex 
interactions between different domains of the Earth sys- 
tem, which may themselves be viewed as complex dynami- 
cal systems, e.g., the atmosphere, hydrosphere, cryosphere 
and biosphere, remains a great challenge for modern sci- 
The urge to make progress in this field is par- 



25 



ence 

ticularly pressing as substantial and mutually interact- 
ing components of the Earth system (tipping elements), 
such as the Indian Monsoon and ENSO, may soon pass 
a bifurcation point (tipping point) due to global climate 
change and consequently experience abrupt and possibly 
irreversible transitions in their dynamics and function |26] . 
Mapping the complex interdependency structure of sub- 
systems, components or processes of the Earth system to 
a network of interacting networks provides a natural, sim- 
plified and condensed mathematical representation. This 
structure could in turn be harnessed for generating new 
insights by graph-theoretical analysis and, hence, foster- 
ing an improved understanding of the Earth system's vul- 
nerability to perturbations like anthropogenic emissions 



of greenhouse gases [27]. As a first step in this direction, 
in this work we develop a novel approach termed coupled 
climate subnetwork analysis for representing and study- 
ing the statistical relationships between several fields of 
climatological observables and apply it to investigate the 
atmosphere's vertical dynamical structure and general cir- 
culation 28,29 from a network perspective. 

So far research on interacting networks has focussed 
on global properties of these systems, e.g., percolation 
thresholds |12| . However, when dealing with networks of 
networks it is of major interest to determine in detail the 
importance and role of single vertices for the interaction 
or communication between different subnetworks as well 
as to characterise their mutual interaction topology, for 
example for studying the vulnerability of coupled and in- 
terdependent networked systems to random perturbations 
or targeted attacks. Here we propose a general and novel 
framework that allows to quantitatively investigate the 
interaction structure of networks of networks on different 
topological scales. We derive graph-theoretical measures 
that allow to answer questions like: Does a vertex have a 
large direct influence on and/or is it an efficient transmit- 
ter of information to a specific subnetwork? Which amount 
of control does a vertex have on the interaction between 
two subnetworks? Are two subnetworks topologically well 
separated or tightly intertwined and is their interaction 
structure well organised or random? The methodology in- 
troduced in this work therefore creates a new setting for 
a detailed graph-theoretical assessment of the functional 
roles of vertices within complex networks of networks, e.g., 
integration or segregation of informatio n in co rticocortical 
networks of cats or macaque monkeys |30||31| . 

After introducing measures and related theoretical con- 
siderations for characterising the topology of interacting 
networks in Sect. [2] we employ this framework in Sect. [3] 
for studying the atmosphere's vertical dynamical struc- 
ture within the context of climate network analysis. Con- 
clusions are drawn in Sect. |U 



2 Theory: The topology of interacting 
networks 

Consider a network G = (V, E), where V — {1, .., N} is a 
set of vertices or elements and E a set of edges or inter- 
actions with N = \V\. As we wish to study a network of 
interacting subnetworks, we consider a decomposition of 
the vertex set V into disjoint sets Vi such that UiVi = V 
and ViHVj = 0,Vi j, where the number of vertices in 
subset Vi is Ni = \Vi\. Similarly, the edge set E is de- 
composed into sets En containing edges between vertices 
inside Vi and sets Eij of edges connecting vertices from Vi 
and Vj, i.e., Uy-Ey = E with E i3 ^E kl = 0, V(i,j) + (fc, I). 
In other words, the mutual interactions between subnet- 
works Gi — (Vi,Eu) (Gi is the induced subgraph of Vi) 
are described by the edge sets E+j for i ^ j. For simplicity 
we restrict ourselves to undirected and unweighted sim- 
ple graphs, since the generalisation of the concepts and 
measures introduced below is straightforward. This type 
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of network is conveniently represented by the symmetric 
adjacency matrix A pq with A pq — 1 if {p, q} £ E and 
A 



i'<i 



otherwise. In this work, indices i,j,k,l always 
denote subnetworks while v,w,p,q designate single ver- 
tices. 



In the following we define several local as well as global 
network measures to quantify and investigate the interac- 
tion topology of networks of networks on different topo- 
logical scales. Here, local network measures assign a real 
number to a vertex v £ V, in relation to (a generally dif- 
ferent) subnetwork Gj, or to any other vertex w £ V de- 
pending on two subnetworks Gi,Gj. They are inspired 
by the "trinity" of classical and frequently used central- 
ity measures degree, closeness and betweenness [2|[3|[32|, 
as well as the local clustering coefficient [33] , and accord- 
ingly quantify direct influence on Gj (cross-degree, Eq. 
(fill), local organisation of interdependency with Gj (lo- 
cal cross-clustering, Eq. |3|), efficiency of interaction with 
Gj (cross-closeness, Eq. T5|)) and the control over com- 



munication between d and Gj (cross-betweenness, Eq. 
(fij)), respectively. The global network measures we in- 
troduce assign a real number to a pair of subnetworks 
(Gi, Gj). They are derived from the well established mea- 
sures edge density, global clustering coefficient and average 
path length 33 . These global network measures quantify 



various overall aspects of the interaction between two sub- 
networks such as its degree of organisation (glob al cros s- 
clustering coefficient and cross-transitivity, Eqs. ( 10|11 1), 
efficiency and speedof information transfer (cross-average 
path length, Eq. (12)) or their mutual interconnectivity 
(cross-edge density, Eq. (|9l) 



Along these lines, similar generalisations of other lo- 
cal and global network properties (like local random walk 
betweenness or global efficiency among many others [l]- 
[3]) may be derived to quantify additional nuances of the 
topology of interacting networks. Adaptations for directed 
and edge- or vertex- weighted networks [34] are also straight- 
forward to deduce. Similarly it is possible on the basis 
of our proposed framework to design measures to take 
into account different qualities or functions of vertices and 
edges within or between subnetworks, e.g., the additional 
constraint that the functioning of vertices v £ Vi depends 
on the functioning of vertices v' £ Vj studied by Buldyrev 
et al. 12 . The selection of measures presented here was 



chosen to be as concise as possible, while at the same time 
representing all classes of commonly used network quan- 
tifiers. 



In the following definitions, we always assume v £ Vi if 
not indicated otherwise. The formula explicitely account 
for the general case i ^ j but can be, nevertheless, eas- 
ily modified to suit the special case i — j. Furthermore, 
the term cross generally relates to the interaction between 
subnetworks Gi,Gj, whereas internal refers to the struc- 
ture within a single subnetwork. 



2.1 Local measures 

Cross-degree centrality k l J gives the number of neigh- 
bours of the vertex v within subnetwork Gj , 



u'i-i — hi — SP A 

"'v v ~ / t -"-vqi 



£ Vi. 



(1) 



This measure thus captures the importance of v for the 
interaction or communication between two subnetworks 
Gi , Gj in terms of the number of direct connections it 
projects between Gi and Gj. For brevity, we will in the 
following suppress the somewhat redundant index i when- 
ever possible, e.g., write k{ instead of k l j . 

The standard degree centrality k v considering the full 
network G can be obtained by summing up the contribu- 
tions from all subnetworks: 



k v — k I . 



(2) 



Local cross-clustering coefficient C l J estimates the prob- 
ability that two randomly drawn neighbours of v from 
subnetwork Gj are also neighbours (Fig. [2]), 



rv — pi — 



ky{kv l) 



AAA 



(3) 



p^qEVj 



For all vertices v* with k J v * £ {0, 1} we set C J V » — 0. C{, 
quantifies the tendency of vertices to form clusters span- 
ning two subnetworks and therefore contains important 
information on the interaction structure between them. 
Assuming no correlations between the occurrence of edges 
within Gj and between Gi and Gj (described by the sets 
Ejj and Eij, respectively) the expectation value for C° v is 
the internal edge density of subnetwork Gj, i.e., (C^) = 
Pj,Vv, where pj = 2\Ejj\/(Nj(Nj - 1)). This lack of cor- 
relations would arise if edges in Eij were distributed ran- 
domly and independently between subnetworks Vi and Vj . 
Hence, C 3 V S> Pj or C 3 V <C Pj indicates significant (anti-) 
correlations in the connectivity between both subnetworks 
pointing at design principles or details of growth processes 
depending on the specific application (e.g., see Sect. 3.3). 



Compared to the other local measures introduced in 
this section, Cl has a less direct relationship with the 
standard local clustering coefficient C v [2], since the con- 
tributions from triangles containing vertices from three 
different subnetworks also have to be taken into account 
(second summand in Eq. (|4|) 



(4) 



fc^-i) fs*^ 1)c - 



^ ] ^ ] A V pApqAq V I , 



where again this holds only if k v > 1. 
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Fig. 2. The local cross-clustering coefficient C l J is the 
probability that two randomly drawn neighbours of vertex 
v from subnetwork Gj are neighbours themselves, where 
v belongs to subnetwork Gi. 



Cross-closeness centrality c l J measures the topological 
closeness of v to subnetwork Gj along shortest paths, 



E 



q<£Vj dvq 



(5) 



where d vq is the shortest path length between vertices v 
and q. If no path exists, i.e., the vertices are not mutually 
reachable, d vq = N — 1 is set as an upper bound, since this 
is the longest possible path length in the considered net- 
work. It is important to note that for generality, we do not 
restrict paths to subnetworks Gi , Gj or a particular order 
of vertices as in earlier works [35] , which might however 
be appropriate depending on the specific application. In 
contrast, the shortest paths analysed here and below may 
contain any vertices w G V in any order depending on the 
topology of the full graph G. 

Cross-closeness therefore quantifies the efficiency of in- 
teraction between a particular vertex and a specific sub- 
network. A vertex with high cross-closeness is likely to be 
important for the fast exchange of information with a cer- 
tain subnetwork, even if it does not have a large number of 
direct neighbours in that subnetwork or a high closeness 
centrality with respect to the whole network. 

The standard closeness centrality c v [32] can be ob- 
tained from the cl as 



N -1 



-1 ■ 



(6) 



Cross-betweenness centrality For any vertex w € V, 
cross-betweenness centrality b 1 ^ indicates its role for medi- 
ating or controlling interactions/communication between 
two subnetworks G, and Gi 



a,, 



^ a. 

p€Vi ,g£ Vj \p,q^w 

= b ji 



(7) 



where o~ pq gives the total number of shortest paths from p 
to q and cr pq (w) counts the number of shortest paths be- 
tween p and q that include w [321136]. While a vertex with 



high cross-degree (relational hub) may function as a robust 
transmitter between two subnetworks, its functional re- 
dundancy for communication is evaluated by cross-between- 
ness centrality. E.g., a relational hub with high cross- 
betweenness is more vulnerable to attack or failure with 
respect to the interaction of two subnetworks than another 
one with low cross-betweenness and, hence, high redun- 
dancy. 

The standard betweenness centrality b w 32 of the full 
network is obtained by summing up the contributions from 
all pairs of subnetworks, 



Ee 



(8) 



In this sense, cross-betweenness centrality can be seen as a 
decomposition of the standard betweenness centrality with 
respect to a certain partition of the full network. Cross- 
betweenness is a generalisation of the measure Q2 defined 
by Flom et al. that is restricted to networks consisting of 
only two different sorts of vertices or two subnetworks in 
our terminology |35}[36| . 



2.2 Global measures 

Cross-edge density pij measures the density of connec- 
tions between distinct subnetworks Gi and G,, 



\E, 



Pij 



NiNj 



(9) 



= Pji- 

Two subnetworks can be considered to be well separated 
topologically if their internal edge densities pi , pj are clearly 
larger than their cross-edge density, i.e., p^ <C pi,Pj- More 
generally, subnetworks form communities of the full net- 
work [13] if this inequality holds for all pairs of subnet- 
works. In. other words, in this situation the partition of 
the full network G induced by the subnetworks Gi would 
give rise to a high modularity 37 . It should be stressed, 
however, that there exists a multitude of other definitions 
of communities [13] . Our general framework does not re- 
quire the chosen partition to be consistent with any such 
definition as long as it allows the resulting cross-network 
measures to be readily interpretable. 



Global cross-clustering coefficient Cij is an estimate of 
the probability of vertices from subnetwork Gi to have 
mutually connected neighbours within subnetwork Gj, 



1 



(10) 



E 

veVi,k J v >i 



Ysp^qEVj A V pA pq A qv 

^2p^ q eVj A vp A vq 



J.F. Donges et al.: Investigating the topology of interacting networks 



5 



It is important to note that in contrast to cross-edge den- 
sity and cross-average path length the cross-clustering co- 
efficient is not a symmetric property of two subnetworks, 
i.e., Cij 7^ Cji. As was shown above, we expect Cy = O(pj) 
if the interaction structure of subnetworks Gi and Gj is 
random, i.e., cross edges are distributed randomly and in- 
dependently between the two subnetworks. In contrast, 
an organised interdependence would induce Cy ^> pj or 
< pj. 












E J 









Cross-transitivity 7y is the probability that two vertices 
in subnetwork Gj are connected if they have a common 
neighbour in subnetwork Gi, 



(11) 



Here also 7y ^ Tji- As the global cross-clustering coeffi- 
cient Cij , Tij is a measure of the degree of organisation of 
the subnetworks' topological interdependence. Analogous 
to the standard version of transitivity [2], Tij tends to 
weigh contributions of vertices in Gj with low cross-degree 
k J v less heavily than Cij. From another point of view, in 
the average of Eq. (Ill edges have equal weights, while Eq. 



( 10 ) has the same property with respect to vertices. This 



peculiarity has to be born in mind when interpreting the 
values of Tij and Cy for general complex networks of net- 
works. Given a fully random interconnectivity structure, 
the expectation value for cross-transitivity is (Tj) = Pj. 



Cross-average path length Cij measures the topological 
closeness of two subnetworks Gi and Gj, 

1,J veVi,qeVj 



= C ■■ 

where My is the number of pairs (v € Vi,q G Vj) which 
are not mutually reachable. Hence, Cij measures the av- 
erage length of existing shortest paths between subnet- 
works Gi and Gj. Cij can be interpreted as a measure 
of the efficiency of interaction between two subnetworks. 
Subnetworks with low Cij are closely interwoven and are 
likely to show a high degree of functional interdependence, 
while those with high dj are topologically more separated 
and likely to be dynamically and functionally more inde- 
pendent of each other. If My = 0, £y is related to the 
cross-closeness centralities c l J via 



C, 



(13) 



3 Application: Analysing the vertical 
dynamical structure of the Earth's 
atmosphere 

Climatologists are interested in studying the correlation 
structure of climatological field variables in order to find 



Fig. 3. Illustration of a coupled climate subnetwork as it 
is constructed in this work, where V\ denotes the set of ver- 
tices in the near ground subnetwork and Vi that of another 
isobaric surface higher up in the atmosphere. En, En are 
sets of internal edges of the two isobaric surfaces or subnet- 
works describing the statistical relationships within each 
isobaric surface, while En contains information on their 
mutual statistical interdependencies. 



spatial as well as temporal patterns accounting for a large 
fraction of the field's variance, commonly relying on estab- 
lished linear methods such as principal component analy 
sis (termed empirical orthogonal functional analysis in cli- 
matological parlance) or singular spectrum analysis 
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|41| . Recently, complex climate networks have been suc- 
cessfully used to analyse single fields of climate observ- 
ables combining methods from nonlinear time series anal- 
ysis and network theory 



17-24 



The formerly mentioned classical approach has been 
generalized to find coupled patterns in climate data in- 
vestigating the correlation structure between fields of dif- 
ferent climatological observables using techniques such as 
canonical correlation analysis or singular value decomposi- 
tion of the covariance matrix |42| . In analogy to the study 
of coupled patterns, here for the first time we expand the 
climate network approach to analyse the dynamical inter- 
relationships between two different climatological fields by 
constructing coupled climate subnetworks and investigat- 
ing them within the interacting networks framework in- 
troduced above. In other words, the correlation structure 
within and between two sets of discrete field observables 
Xi (t) , Yj (t) is mapped to a complex coupled climate sub- 
network (Fig. [3]), where i and j are spatial indices and t 
denotes time. 

In order to justify treating a coupled climate subnet- 
work as a network of networks, an ab initio physical sep- 
aration of the climatological fields is necessary regarding 
processes responsible for internal coupling within a sin- 
gle field and those mediating interactions between both 
fields, as is elaborated in the following. The Earth's quasi- 
spherical shape and almost homogeneous mass distribu- 
tion result in a hydrostatic equilibrium in first-order ap- 
proximation, implying a stable isobaric quasi-horizontal 
stratification and therefore a strong buoyancy constraint 
[43) . Local heating of the Earth's surface and atmosphere 
due to temporally and spatially varying solar radiation 
induces minor disturbances of the system, which give rise 
to weather variability and are propagated by advection, 
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diffusion (dominantly turbulent in the homosphere) and 
convection processes. Convection processes lead to verti- 
cal movement resulting in pressure gradients which are 
balanced by quasi-horizontal geostrophic winds along iso- 
bares. An important mechanism is slant convection arising 



as an adjustment of baroclinic instabilities 44 . 

Geopotential height Z(j), <j>, h) is a vertical coordinate 
referenced to the Earth's mean sea level which takes into 
account the variation of gravitational acceleration g(j) 1 <fi, h) 
with latitude z9, longitude 4> and geometrical height h. It 
is defined as 

Z{#,J),h) = - ( g(#,cj>,z)dz, (14) 
9o Jo 

where go denotes the standard gravitational acceleration 
at mean sea level [43] , Z(z9,</>, /i) is approximately equiv- 
alent to geometrical height h within the homosphere, i.e., 
the lower portion of the atmosphere we consider in this 
work. The geopotential height Z{&, </>, P) of a certain pres- 
sure level P is defined as the geopotential height necessary 
to reach the given pressure P. In meteorology and clima- 
tology, the field Z(d, cf>, P) is frequently used as an equiva- 
lent and convenient representation of the three-dimensional 
atmospheric pressure field P(i), <j>, Z). Therefore the dis- 
cretised and vertically resolved geopotential height field 
Zl (t) sampled at predefined points v on isobaric surfaces i 
at pressure P t captures the dynamics of both the geostrophic 
wind field as well as convection processes and, hence, re- 
flects global weather and climate dynamics to a good ap- 
proximation [43] . Given the distinct physical processes be- 
hind vertical and quasi-horizontal atmospheric dynamics 
described above, it is hence feasible to apply the interact- 
ing networks approach, treating the induced subgraphs of 
vertices lying on the same isobaric surface i as distinct 
subnetworks Gj. In this paper we specifically focus on the 
interaction structure between near ground and upper level 
atmospheric dynamics, which is particularly interesting as 
a large portion of the solar forcing driving atmospheric dy- 
namics takes place on the Earth's surface (Fig. [3]). 

To better understand our coupled climate subnetwork 
approach, one should be aware of the strong analogy ex- 
isting between the concepts of climate networks and func- 
tional brain networks studied in neuroscience and medicine 
6,7, 30 , 31 , 45 . Both types of networks describe statistical 
similarity relationships between spatially embedded time 
series using the same methods of time series analysis, but 
relying on distinct data sources, i.e., the climate system 
and the mammalian brain. Now there exist two fundamen- 
tally different types of networks: (i) Structural (or anatom- 
ical) networks, on the one hand, reflect the topological 
structure of existing ties between objects (e.g., computers, 
neurons, columns of neurons), referring to either physi- 
cal connections (e.g., internet, power grids, neuronal net- 
works) or abstract relations (e.g., world wide web, social 
networks, citation networks), (ii) Functional networks, on 
the other hand, including functional brain networks and 
complex climate networks, are extracted from an under- 
lying system by detecting and assessing similarities in the 
dynamical behaviour of its components. In other words, 



structural networks represent a priori knowledge on a sys- 
tem's internal structure on a certain level of abstraction, 
whereas functional networks are inferred solely from the 
measured or simulated dynamics of subsystems, usually 
without including any additional information. Hence, in 
contrast to structural networks, the resulting topological 
interconnections in functional networks do not directly al- 
low to draw conclusions on a causal interrelationship be- 
tween vertices. When constructed and compared for the 
same system, e.g. a neuronal network, both types of net- 
works will usually not be identical. This implies that spe- 
cial emphasis has to be put on physical arguments when 
interpreting topological features of climate networks. A 
further noteworthy duality exists between spatial similar- 
ity networks based on fields of time series such as the cli- 
mate networks discussed above, and temporal similarity 
networks, e.g., recurrence networks 
a single time series. 



46 48 



representing 



3.1 Data 

To construct coupled climate subnetworks capturing longer- 
term dynamics of the geostrophic wind field as well as 
large-scale convection processes, we rely on the global month- 
ly averaged and vertically resolved atmospheric geopoten- 
tial height field covering the troposphere and the lower 
stratosphere. We use Reanalysis 1 data provided by the 
National Center for Environmental Prediction/National 
Center for Atmospheric Research (NCEP/NCAR) [49). 
For each of the 17 isobaric surfaces P l the NCEP/NCAR 
data is given on an equally spaced spherical grid with a lat- 
itudinal and longitudinal resolution of 2.5° x 2.5°, resulting 
in 10, 266 grid points. Using this type of grid for network 
construction would induce biases in the statistical prop- 



Table 1. Pressure and associated mean geopotential 
height Zi (Eq. (15)) for each isobaric surface i in the 
NCEP/NCAR Reanalysis 1 reconstruction of the geopo- 
tential height field. 



/ 


Pi [mbar] 


Zi [km 


1 


1000 


0.1 


2 


925 


0.8 


3 


850 


1.5 


4 


700 


3.0 


5 


600 


4.3 


6 


500 


5.7 


7 


400 


7.3 


8 


300 


9.3 


9 


250 


10.6 


10 


200 


12.0 


11 


150 


13.8 


12 


100 


16.3 


13 


70 


18.5 


14 


50 


20.5 


15 


30 


23.8 


16 


20 


26.4 


17 


10 


30.9 
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erties of climate networks, since the area covered by each 
grid point is not uniform but decreases towards the poles 
like the cosine of latitude [23][34] . To avoid these effects, 
we choose to project the data to an icosahedral grid [50] 
of N, t = 2, 562 time series Z*(t) with v G {1, . . . , iV;} for 
each isobaric surface i, respectively, at pressure Pi using 
the conservative interpolation scheme described in [51] . 
Each isobaric surface i may be associated to an average 
geopotential height 



W*)>„« 



(15) 



by averaging over time t and all grid points v contained 
within this level (Tab. [TJ. The time series Z l v {t) contain 
734 monthly averaged data points from January 1948 until 
February 2009. 

Relying on the icosahedral grid, which is used by the 
German Weather Service for its operational global weather 
forecast model Global Modell Extended (GME) and re- 
ferred to as "triangular grid" 
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guarantees nearly uni- 
form grid cell areas within each isobaric surface i. The 
differences in grid cell area between different isobaric sur- 
faces due to their varying distance from the Earth's sur- 
face are negligible since the maximum vertical separation 
of isobaric surfaces (« 30 km, see Tab. [I} is much smaller 
than the Earth's mean radius a ~ 6, 370 km. 

As a final step of preprocessing the data, we calculate 
climatological anomaly time series Z % v (f ) by phase averag- 
ing (see [l8]) to remove the leading order effect of the an- 
nual cyclefrom the geopotential height time series Z l v {t). 
This helps to avoid spurious correlations solely due to the 
solar forcing common to all time series. 



{z p (t),z q (t)} i 



is then given by 



p*ii 



(Z p (t) -Ai p ) {^q{t) -f J 'l)) t 



(16) 



where [i pq and <r Ptq are the mean and standard devia- 
tion of the Z p ^ q (t) , respectively. Edges {p, q} are intro- 
duced into the coupled climate subnetwork if the abso- 
lute value of Pearson correlation \P1'^\ at zero lag be- 
tween both time series Z p {t) and Z q {t) exceeds a cer- 
tain threshold value < T < 1. The threshold should 
be carefully chosen to ensure that only statistically signif- 
icant and reasonably strong correlations are included in 
the network |18| . Because an optimal value of T is neither 
easy to define nor to determine (it would even be possible 
to determine a threshold T pq for each pair of time series 
p, q separately), we w ill present results for various thresh- 
olds below (Sect. 3.3). Furthermore it would be feasible to 



avoid the choice of a threshold altogether by including a 
function of the correlation measure w p '-> = f{P p ^) into the 
network analysis as edge weights w p ^ . However, we do not 
follow this research avenue here as it introduces as a new 
complication the choice of a meaningful transfer function 
/ for mapping correlations to edge weights depending on 
the interpretation of a specific network measure. In con- 
trast, this work focusses on the topology of statistical in- 
terrelationships between different isobaric surfaces on the 
network level. The consideration of edge weights will be 
an interesting subject of future work. 

Finally, the coupled climate subnetworks are completely 
described by the adjacency matrices 



3.2 Network construction 

We now construct a sequence of pairwise coupled climate 
subnetworks based on statistical interrelationships within 
the three-dimensional geopotential height anomaly field 
by using the linear Pearson correlation at zero lag. As 
was shown 



18 



linear correlation measures are suf- 
ficient for a first overview study like the present work, 
because they capture the great majority of statistical in- 
terrelationships within fields of smooth (non-intermittent) 
climatological variables like temperature or geopotential 
height. The so obtained networks describe both the in- 
trinsic structure of a single isobaric surface as well as 
the interaction structure between two isobaric surfaces, 
in other words comprising horizontal as well as vertical 
interdependencies of the spatially embedded time series. 
First the Ni,Nj = 2,562 time series of two isobaric sur- 
faces i,j are relabeled using indices p,q G {1,...,A}, 
N = Ni+Nj, Z l v (t),Zl(t) ->■ Z p (t),Z q (t), where the fully 
coupled climate subnetworks contain N = 5, 124 vertices. 
The anomaly time series Z p , q (t) are identified with vertices 
p, q of the coupled climate subnetwork and subsequently 
collected in sets Vi,Vj based on their native isobaric sur- 
face. The classic Pearson correlation between any pair of 
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where ©(•) is the Heaviside function and Kronecker's delta 
S pq indicates that artificial self-loops are not considered. 
The subnetworks Gi , Gj representing the correlation struc- 
ture within each isobaric surface i , j are the induced sub- 
graphs of the sets Vi,Vj embedded within the coupled 
climate subnetwork G lJ described by the adjacency ma- 
trix A p ^. In the following, local as well as global (cross- 
network measures m l J,mij will be calculated from the 
coupled climate subnetwork G l ' J consisting of two isobaric 
subnetworks Gi and Gj and their interaction structure 
i.e., G 1 ' 3 = (Vi U Vj,Eu U Ejj U E l3 ). 
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3.3 Results 

Global measures The available data describes the dy- 
namics of geopotential height in the lower homosphere, en- 
compassing the troposphere and lower stratosphere, where 
most atmospheric dynamical processes relevant for the 
Earth's climate system are concentrated to 43 . In the 



following we will investigate which aspects of atmospheric 
dynamics can be revealed by analysing the sequence of 



coupled climate subnetworks G 



l.i 



1,...,17. 
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A B 




Fig. 4. (A) Internal edge density pi of all subnetworks on isobaric surfaces i, (B) cross-edge density pu between the 
surface level 1 and all other isobaric surfaces i of average height Zi, (C) the ratio pi/pu and (D) the ratio pi/pu in 
the coupled geopotential height climate network for various thresholds T. Note the inversions of cross-edge density at 
Z m 1 — 3, 12, 16 and 26 km becoming increasingly pronounced for decreasing threshold T. 



It is pivotal to be aware that for similarity networks 
constructed from spatio-temporal data, such as the cou- 
pled climate subnetworks studied here, network structure 
is subject to statistical uncertainties [53] - One consequence 
is that average path length and global clustering coeffi- 
cient cannot be considered as useful order parameters for 
network classification 
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Both are biased by spurious 
or missing links due to the network construction algo- 
rithm and local correlations between spatially close ob- 
servations. Particularly, this implies that it is not mean- 
ingful to classify interaction networks as small-world net- 
works [33| , Erdos-Renyi type random networks or grid-like 
regular networks, as is common practice for other types 
of networks [2} [3]. S ince the corresponding arguments of 
Bialonski et aL~|45| also hold for the modified clustering 
coefficients, transitivities and average path lengths defined 
in Sect.[2j we in the following do not consider the absolute 
values of these measures. In contrast, here we are solely 
interested in their relative changes conveying information 
on the varying interaction structure within and between 
different isobaric surfaces of the atmosphere. 

Our first and fundamental observation is that the cross- 
edge density pu (Fig. |4j3) between the near-surface and 
all other isobaric surfaces i is always smaller than and 
well separated from the internal edge density of the up- 
per isobaric subnetwork pi (Fig. |4|\) for all considered 



thresholds T. Moreover, the ratio of the edge densities 
Pi/pu and, hence, the physical separation of the underly- 
ing dynamics is increasing with height Zi (Fig. [J)). Ap- 
proximately the same holds for the internal edge density 
pi of the near-surface subnetwork which is considerably 
larger than pu for most Zi and T (Fig.|4]C). This observa- 
tion reflects topologically the dynamical separation of at- 
mospheric processes within and between isobaric surfaces 
i, respectively, that led to identifying them with subnet- 
works in the first place (see the introduction to Sec. [3]). 
Nevertheless, p\ can be slightly smaller than the cross- 
edge density pu, particularly close to the two pronounced 
minima of the ratio p\ / pu (the overall minimum value is 
mm^r pi/pu ~ 0.7, see Fig. (4b). This finding, however, 
does not challenge the applicability of our method as our 
framework does not require that subnetworks constitute 
communities of the full network (Sec. 2.2). 



The cross-edge density pu displays prominent extrema 
with varying height Zi of the isobaric surface i, which be- 
come more pronounced for decreasing threshold T (Fig.[4|3) 
Two maxima of pu are located at 1—3 km and 16 km, while 
a minimum is found at 12 km. A much more weakly devel- 
oped inversion of cross-edge density pu occurs at 26 km, 
but it is only visible for small T. These findings indi- 
cate that correlations of large-scale quasi-geostrophic wind 
dynamics are significantly increased between the near- 
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Fig. 5. The cross-average path length Cu between the 
surface level 1 and all other isobaric surfaces i of average 
height Zi in the coupled geopotential height climate net- 
work for different thresholds T. In the lower stratosphere 
and for larger thresholds, Cu is not defined due to vanish- 
ing cross-edge density pu and, hence, the corresponding 
data points are not shown. Cu shows inversions at similar 
height levels as the cross-edge density pu (Fig. [4| , which 
in contrast decrease in sharpness for deceasing threshold 
T. 



ground and higher isobaric surfaces at approx. 1—3 km and 
16 km. The superficially similar inversions in internal edge 
density pi with two maxima at 4 — 6 km, 18 — 24 km and a 
minimum at 7 — 1 1 km have to be carefully distinguished 
from those of cross-edge density pu (Fig. |4]A_). First, the 
geopotential height intervals within which the respective 
extrema are observed for different T do not overlap. Sec- 
ond, on the one hand, recall that in contrast to pu the in- 
ternal edge density pi measures dynamical correlations oc- 
curring within a quasi-horizontal isobaric surface at pres- 
sure Pi. On the other hand, the physical processes act- 
ing within isobaric surfaces (geostrophic wind, planetary 
Rossby waves, gravity waves) and between them (convec- 
tion, turbulent mixing), which are relevant for large scale 
dynamical coupling are distinctively different. 

The cross-average path length Cu between the near 
ground layer and all other isobaric surfaces possesses two 
minima at 3 km and 16 km as well as one maximum at 
If km (Fig. [5]). In contrast to cross-edge density, small 
values of Cu imply tight dynamical relationships between 
two isobaric surfaces, while large values indicate a weaker 
coupling. Hence, cross-average path length consistently 
behaves complementarily to cross-edge density pu and in- 
ternal edge-density pi, as shortest paths between two dif- 
ferent subnetworks generally contain edges from both sub- 
networks implying their average length to decrease with 
increasing pu and pi. Interestingly, in contrast to the be- 
haviour of pu and pi, the extrema in Cu are rendered 
increasingly pronounced for increasing threshold T . Cross- 
average path length remains sensitive to variations in the 
topological closeness of two isobaric surfaces even as more 
and more edges (p, q) of weaker correlation strength |P pg | 
are removed from the coupled climate subnetworks. 



While the cross-network measures discussed so far are 
symmetric with respect to exchanging the involved sub- 
networks (see Sect. [2]), the two clustering measures global 
cross-clustering coefficient Cij and cross-transitivity Tij 
are intrinsically directional. Hence, as for the horizontally 
resolved cross-measures to be treated below, we are able 
to distinguish Cu and Tu pointing "upward" from the 
near ground to higher isobaric surfaces from their coun- 
terparts Cu and Tu projecting "downward" . Both Cu and 
Tu consistently uncover that the probability of vertices 
within the near ground isobaric surface to have connected 
neighbours in higher isobaric surfaces reaches local max- 
ima between 3 — 6 km and 14 — 16 km, whereas a local 
minimum is assumed at 11 km (Figs. [6]^ and C). For both 
measures, all three inversions are only observed for low 
thresholds T. It is particularly interesting to compare the 
measured values of Cu and Tu with those expected for a 
fully random c onne ctivity structure between isobaric sur- 
faces (see Sect. 2.2). The "upward" pointing global cross- 
clustering coefficients Cu are considerably larger than the 
expectations pi (indicated by dashed lines in Fig. [6]A_) be- 
low 15 km, and markedly smaller for average geopotential 
heights above 20 km. Similarly, the "upward" projecting 
cross-transitivity Tu is significantly larger than the ex- 
pected values pi for the fully random null model for all 
height levels (Fig. [6p) . 

Compared to their counterparts the behaviour of the 
"downward" projecting measures Cu and Tn is less consis- 
tent. The global cross-clustering coefficient Cu possesses 
two extrema at 3 — 6 km and 16 — 21km, whereas the 
cross-transitivity Tu takes local minima between 2 — 3 km 
and 14 — 16 km and local maxima between 12 — 14 km and 
at 24 km of geopotential height (Figs. [6j3 and D). Note 
that "downward" pointing cross-transitivity Tn behaves 
complementarily to its "upward" projecting counterpart 
Tu- Both clustering measures are significantly larger than 
their values p\ expected for a fully random connectivity 
structure for nearly all height levels, only above 20 km the 
global cross-clustering coefficient Cu is consistent with the 
expectation value for larger thresholds T. 

In summary, the clustering measures exhibit that the 
interaction topology between the near ground and higher 
isobaric surfaces is not consistent with a fully random null 
model, except for the lower stratosphere above 20 km. The 
same holds when comparing the observed values of the 
clustering measures to those expected from a more so- 
phisticated null-model with fixed cross-degree sequences 
but otherwise random interaction topology (results not 
shown here) . These findings highlight that the coupled cli- 
mate subnetworks considered here have a nontrivial inter- 
action topology, which is consistent with known features of 
the atmosphere's vertical dynamical structure and strati- 
fication and may contain additional information on previ- 
ously unknow n feature s of t he atmosphere's general circu- 
lation 28,29 (see Sect. 3.4). The large values of clustering 
measures observed for most height levels can be partly 
explained by spatial correlations between closely neigh- 
boured time series stemming from the continuity of the 
geopotential height field in conjunction with the intrinsic 
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Fig. 6. Global cross-clustering coefficients (A) Cu calculated "upward" from the surface level 1 to all other isobaric 
surfaces i and (B) Cu defined "downward" (continuous lines), as well as cross-transitivities (C) Tu ("upward") and 
(D) Tn ("downward"). For comparison with the global cross-clustering coefficients (Cj,-) expected for a fully random 
connectivity between isobaric surfaces, (A) and (B) also feature the expectation values (Cu) — pi and (Cu) = pi, 
respectively (dashed lines). 
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54 . Directly comparing Cu 
i, respectively, furthermore 
clearly reveals the bias in the global cross-clustering coeffi- 
cient, which leads to generally lower values for increasing 
height through weighing more strongly the contribution 
of the increasing number of vertices of low cross-degree 
induced b y th e decreasing trend in cross-edge density pu 
(see Sect. 2.2). Consequently, as for the standard versions 
of global clustering and transitivity [2], special care has 
to be taken when interpreting absolute values of global 
cross-clustering and cross-transitivity. The suggested best 
practice is to always consider the two measures simultane- 
ously and to draw conclusions only from qualitative fea- 
tures exhibited by both of them. 



Local measures For visualising the inherently three-dimen- 
sional fields of local cross-network measures (one of the 
subnetwork indices i,j indices being fixed as in our ap- 
plication) m l v J = rn£,g ^ , where i? and <j> denote latitude 
and longitude, we choose to focus on their variation with 
height and latitude. This is most appropriate as our aim 
is to study from a complex network perspective aspects of 



transitivity of the Pearson correl ation coefficient P pq used the atmosphere's general circulation 28|29 , the dominant 
for network construction 
and Tn as well as Cu and 



forcing of which is the latitudinal variation of radiative 
solar forcing 43 . Hence, in the following we will consider 



zonal averages 



(18) 



along circles of constant latitude. A detailed study and 
interpretation of the latitudinally and longitudinally re- 
solved fields of cross-network measures will be the subject 
of future work. Like the scalar measures of cross-clustering 
discussed above, most local cross-network measures are 
non-symmetric with respect to interchanging subnetworks 
and, hence, intrinsically directional. Therefore we distin- 
guish "upward" cross-degree centrality k u (fl) from "down- 
ward" cross-degree centrality fc ll (i?), and "upward" cross- 
closeness centrality c 11 ^) from "downward" cross-close- 
ness centrality c ll ($). For brevity here we present results 
for representative thresholds of T = 0.4,0.5,0.6 only. 

Cross-degree and cross-closeness centrality show a sim- 
ilar structure in both directions (Figs. [7j [8j [^V-D). Both 
measures are generally increased in the tropics and po- 
lar regions, whereas they take smaller values in the mid- 
latitudes. We observe a pronounced asymmetry between 
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Fig. 7. Zonally averaged cross-degree centralities (A) k 11 ^) pointing "upward" from the near ground level 1 to 
all other isobaric surfaces i and (B) fc ll (t?) projecting "downward", zonally averaged cross-closeness centralities (C) 
c ll ($) pointing "upward" and (D) c ll (t7) projecting "downward", (E) near ground and (F) b\ l (-d) upper level 

component of zonally averaged cross-betweenness centrality for a threshold of T = 0.4. Panel (B) can be interpreted 
to show the number of cross-edges connecting a certain volume element with the whole near ground isobaric surface, 
averaged along bands of approximately equal latitude (approximately because of the geodesic grid). 



both hemispheres as the cross-degree and cross-closeness 
centralities in the northern polar regions are significantly 
larger than those above Antarctica and the surrounding 
Southern Ocean in the southern hemisphere. Furthermore, 
a structure of increased "downward" cross-degree fc ll ($) 
and -closeness centrality d l (d) appears above the north- 
ern mid-latitudes but not above those of the southern 
hemisphere (Figs. Pf) M [9b and D). 



Additionally considering the dependence of both local 
measures on geopotential height Zi, the "upward" point- 
ing cross-degree k 1 " 1 ^) and -closeness centralities c ll (d) 
possess two distinct tropical maxima centred around 1 km 
and 16 km, as well as two northern polar maxima at 3 km 
and 16 km (Figs. ff\ § |K and C). While k u (d) decreases 
monotonously with height in the mid-latitudes and above 
the polar regions of the southern hemisphere, c ll {d) main- 
tains its bimodal structure there. The "downward" pro- 
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Fig. 8. As Fig. with threshold T = 0.5. 



jecting cross-degree k 11 ($) and -closeness centralities c l1 ($) 
display two pronounced tropical maxima at 3 km and 16 km, 
and a maximum centred around 14 km in the northern 
mid-latitudes (Figs. 0JU and D )- Furthermore, c ll (tf) 
reveals a maximum oitopological closeness to the near 
ground isobaric surface at 1 km above the northern polar 
regions, while k ll ($) decreases monotonically with height 
there. 



It is worth noting that the observed cxtrema in cross- 
degree are directly related to those in cross-edge density 
via Eq. Q and extrema in cross-closeness have an equiva- 
lent association to those of cross-average path length (Eq. 
(13)). Moreover, we point out that plots of cross-degree 
centrality like those presented in this work may be used 



to draw conclusions on the main sources and destinations 
of cross-edges without relying on full three-dimensional vi- 
sualisations of the coupled network structure. For exam- 
ple, consider the region of increased "downward" cross- 
degree fc ll ($) between 11km and 18 km of geopotential 
height above the northern mid-latitudes (Fig. [7j [8j [9)3) ■ 
It implies that a considerably large number of cross-edges 
connect this region directly to the whole near ground iso- 
baric surface. We learn where exactly on the near ground 
surface those cross-edges originate from by looking at the 
"upward" cross-degree fc ll ($). It measures how many and 
from where on the near ground surface cross-edges con- 
nect to some higher isobaric surface. Now the regions of in- 
creased fc u ($) between 11 km and 18 km above the tropics 
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and northern polar regions imply that many cross-edges 
originating in the tropics and northern polar regions of the 
near ground surface project to isobaric surfaces between 
11 km and 18 km (Fig.[7j[8l[9k). As these are the only ma- 
jor structures in this range of geopotential height, we may 
conclude that a significant number of cross-edges must link 
the tropical and northern polar near ground isobaric sur- 
face with the northern mid-latitudes' upper troposphere 
to lower stratosphere between 11km and 18 km. 

In contrast to the latter two local measures, cross- 
betweenness b 1 ^ with w G V^U Vj ; is symmetric with respect 
to exchanging the involved subnetworks (see Eq. Q), but 
assigns a value to vertices of both subnetworks. Therefore 
we will in the following analyse zonally averaged fields of 



cross- b et weenness 

* W -M«< (19> 

for vertices taken from a specific isobaric subnetwork i. It 
was shown earlier for climate networks constructed from 
surface air temperature data that betweenness centrality 
(Eq. ([8])) yields additional information when compared to 
degree (Eq. Q) and closeness centrality (Eq. §6§) [Tf) . 
Similarly, the near ground and higher isobaric surface com- 
ponents of cross-betweenness centrality bY^d) re- 
veal rich structures which are partially complementary to 
those seen in the zonally averaged fields of degree and 
closeness centrality (Figs. [7j [8j (9^5 and F). Both fields 
highlight how frequently certain regions on the two iso- 



14 



J.F. Donges et al.: Investigating the topology of interacting networks 



baric surfaces are traversed by shortest paths connect- 
ing the near ground to a higher isobaric surface. Now 
due to our network construction procedure (Sect. 3.2), 



nect ed impose s a substantial constraint on network topol- 
ogy 
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shortest paths correspond to sequences of strongly and 
significantly statistically interrelated pairs of time series 
with a minimum number of intermediate steps. Hence, it 
is conceivable to assume that to a first approximation a 
markedly increased cross-betweenness centrality indicates 
that a region is particularly important for mediating in- 
teractions between two isobaric surfaces, while the con- 
trary is true for regions with significantly decreased cross- 
betweenness [17] . 

The hemispherically asymmetric near ground compo- 
nent of cross-betweenness centrality reveals that 
the northern subpolar regions are exceptionally important 
for mediating interactions between the near ground and 
all considered heights ranging from the lower troposphere 
across the tropopause to the lower stratosphere (Figs. [7] 
[8] [9j3) . The near ground tropics and southern polar re- 
gions appear only to be relevant for coupling the near 
ground to isobaric surfaces in the upper troposphere and 
above. In contrast, the upper surface component of cross- 
betweenness (1?) possesses a more pronounced hemi- 
spherical symmetry (Figs. [7] [8] ^^)- Most notable and 
stable for various thresholds T are the two tongues of in- 
creased b\ l (d) ranging from the near ground tropics to 
the mid-latitudes and subpolar regions in the lower strato- 
sphere. These structures indicate that the latitude at which 
shortest paths connecting near ground and upper isobaric 
surfaces arrive in the upper surface tends to increase to- 
wards the poles for growing geopotential height Zj in both 
hemispheres. One should be aware that the structures de- 
tected in both components of the cross-betweenness field 
in the stratosphere as well as in all other local and global 
measures above 20 km of geopotential height should be 
treated with care. This is because only very few edges ex- 
ist between the near ground and isobaric surfaces in these 
height levels (Fig. ]4j3) and, hence, statistically expected 
false detections or omissions of cross-edges are likely to in- 



duce recognisable changes in the respective measures 45 



To quantitatively assess the effects of small changes in the 
interaction topology between subnetworks, new types of 
significance tests based on network null models need to be 
developed in future work. 

After describing the results of our analysis we would 
like to draw attention to the fact that the observed corre- 
lations in various measures revealed by qualitatively sim- 
ilar structures like the inversions at mostly three differ- 
ent height levels are not necessarily a direct consequence 
of the measure's definitions (see Sect. but point to a 
specific type of network structure. While correlations in 
different measures quantifying distinct aspects of network 
topology need not be present for general network struc- 
tures, they are prevalent in many different types of real- 
world networks and network models 4] . Particularly, these 
correlations are expected to arise in spatially embedded 
functional climate (and brain) networks like the coupled 
climate subnetworks considered in this work, since the in- 
creased probability of spatially close vertices to be con- 



3.4 Climatological interpretation 

We are now in a position to elaborate on the climatological 
implications of our coupled climate subnetwork analysis. 
First, recall that the coupled climate subnetworks were 
constructed from monthly averaged time series of geopo- 
tential height describing the dynamics of the atmosphere's 
quasi-geostrophic wind field on longer than monthly time 
scales. Variability on shorter time scales, e.g., synoptic 
scale weather systems, is included in the averages but does 
not appear explicitly in the time series. Therefore, we can 
indeed expect the coupled climate subnetworks to rep- 
resent the climatological mean state of the atmosphere's 
three-dimensional correlation structure, excluding the di- 
rect effects of such weather phenomena with typical life- 
times of clearly less than one month. 

The cross-network measures discussed in Sect. [32] re- 
veal for a wide range of thresholds T aspects of the at- 
mosphere's stratification and, more importantly, physi- 
cal processes which couple the dynamics on different iso- 
baric surfaces despite the strong buoyancy constraint im- 
posed by vertical stratification. First it should be noted 
that in accordance with physical considerations and obser- 
vations, all cross-network measures consistently indicate 
that most atmospheric dynamics takes place within the 
troposphere with comparatively weak coupling to the su- 
perjacent stratosphere [43]. Steadily increasing (decreas- 
ing) from the near ground to reach a first maximum (min- 
imum) at 1km to 3 km of geopotential height, cross-edge 
density pu and cross- average path length Cu indicate that 
the near ground and higher isobaric surfaces become dy- 
namically more densely interwoven when ascending from 
the planetary boundary layer below approx. 1 km into the 
lower free atmosphere (Figs. [4j3 and [5 ) . Bearing in mind 
the spatial continuity of the geopotential height field, spa- 
tially close vertices are likely to be connected and to share 
common neighbours. This implies that the typical corre- 
lation radius of near ground vertices to higher isobaric 
surfaces and vice versa increases throughout the lower 
troposphere. This observation is consistent with the in- 
fluence of the Earth's surface orography on atmospheric 
flow exponentially decreasing with height in the planetary 
boundary layer via the Ekman effect. Hand in hand with 
the prevalence of turbulence, the formation of long-range 
dynamical couplings is inhibited by this essentially fric- 
tional effect. In turn, within the less turbulent free atmo- 
sphere above approximately 1 km the wind field behaves 
quasi-geostrophically and allows for the long-range prop- 
agation of dynamical influences between the near ground 



and higher altitudes 43 



Within the cross-edge density (Fig. |4p), we observe 
that the first maximum shifts from approximately 4 km 
to 3 km for higher thresholds. This detail can be credited 
to the decreasing height of the planetary boundary layer 
with latitude, the typical zone within which convection 
cells form. The prominent tropical convection processes 
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are known to be more diffusive than those in the mid- 
latitudes and, hence, yield lower correlation values. Thus 
higher thresholds rather account for mid-latitude convec- 
tion phenomena rejecting the relatively low correlated pro- 
cesses in the tropics. The prevalence of the latter, on the 
other hand, is reflected in lower thresholds. 

Above approx. 3 km of geopotential height, pu and 
Cu decrease (increase) again until reaching a local mini- 
mum (maximum) between 11km and 14 km (Figs.|4j3 and 
pj), indicating a dominating effect of a more stable ver- 
tical stratification acting to inhibit dynamical couplings 
between the now vertically more separated isobaric lev- 
els 1 and i. This may be understood considering that 
turbulent vertical mixing is significantly less prevalent in 
the free atmosphere than it is in the planetary bound- 
ary layer and convection is suppressed by baroclinic ad- 
justment forcing the atmosphere to interact horizontally. 
The second local maximum (minimum) of pu and Cu at 
16 km highlights that the cumulative action of tropical 
penetrative convection processes (hot towers) reaching up 
to this height mediates markedly increased dynamical in- 
terrelationships between the n ear g round and isobaric sur- 
faces close to the tropopause [56] . The tropical origin of 
this coupling is more readily seen in the fields of zonally 
averaged cross-degree and -closeness centralities (Fig. [7] 
[8j [9| . When entering the stratosphere, quickly decreasing 
(increasing) pu and Cu indicate that influences reaching 
from the near ground to these heights are strongly inhib- 
ited by the dynamical barrier formed by the temperature 
inversion which marks the boundary between troposphere 
and stratosphere. This conclusion is further supported 
by the random-like interconnectivity structure revealed 
by the "downward" pointing global cross-clustering coeffi- 
cient Cu within the stratosphere (Fig.[6j3). The behaviour 
of internal edge density pt is consistent with the forego- 
ing argumentation in the troposphere, its striking increase 
above the tropopause reveals the spatially uniform dynam- 
ics of the statospheric wind field (Fig. |4]4) [43] . 

The zonally averaged fields of cross-degree, -closeness 
and -betweenness centrality uncover fea tures o f the atmo- 
sphere's general meridional circulation 28 29 . Most evi- 



dent are the Hadley and polar cells which are indicated by 
markedly increased values of both measures in the tropics 
and polar regions (Figs. [7] [8] [9|. Here the generally rising 
motion of air above the equator and subpolar latitudes 
couples surface wind dynamics to the upper troposphere, 
which becomes apparent in both components of zonally 
averaged cross- betweenness (Figs. [7j [§ [9^ and F). The 
surface component of cross-betweenness may also be inter- 
preted to show a signature of the northern hemisphere cir- 
cumpolar vortex around a typical height of approximately 
5 km which is known to induce vertical air motion and, 
hence, vertical dynamical coupling (Figs. [7] [8] [9^3) [57] . 
Supporting this interpretation, a corresponding signature 
is not seen in the southern hemisphere which is consistent 
with the Antarctic ice shield inhibiting the formation of 
a polar vortex there. Cross-degree and -closeness central- 
ity also show that the height of the tropopause decreases 
towards the poles by a slight poleward shift of their max- 



imum values towards lower altitudes, as this dynamical 
barrier strongly inhibits the propagation of signals from 
the near ground to the stratosphere. The Ferrel circulation 
may be involved in forming the comparatively weak dy- 
namical interrelationships between the near ground south- 
ern subtropics and mid-latitudes with isobaric surfaces 
lying in the upper troposphere and lower stratosphere 
(Figs. [7] [8] [9]^ and C). Similarly, the remarkable coupling 
of wind dynamics within the upper troposphere and lower 
stratosphere of the northern mid-latitudes with the near 
ground tropics and northern polar regions might be re- 
lated to the northern Ferre l cel l (Figs.[7][8l[9l3 and D, see 
also the discussion in Sect. 



3.3) 



3.5 Outlook 

The analysis performed in this section serves as an illus- 
tration of the potentials of our network-based approach 
for studying statistical interrelationships between differ- 
ent fields of climatological observables as well as arbitrary 
data with a similar structure. It was shown that the se- 
quence of coupled climate subnetworks constructed from 
three-dimensional geopotential height data contains a lot 
of information on known features of the atmospheric gen- 
eral circulation and stratification 58 . However, partic- 



ularly the measure cross-betweenness centrality is read- 
ily interpretable in the context of coupled climate sub- 
networks (see above), but reveals interesting structures, 
which we cannot obviously relate to known features of the 
atmosphere's dynamical structure. This in turn indicates 
that coupled climate subnetworks in general and cross- 
betweenness centrality in particular have the potential to 
uncover previously unknown features of the atmosphere's 
general circulation and may prove useful in the future to 
shed light on so far open questions of atmospheric dynam- 
specifically those considering the response to cli- 



2<s 
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mate change [26] . While it should be born in mind that the 
results of our study are subject to known deficiencies and 
limitations of the NCEP/NCAR Reanalysis 1 data, the 
variable geopotential height we analyse here is considered 
as one of the most reliable products of this reanalysis, since 
it is strongly determined by observations and, hence, less 
dependent on the particular model used for data assimila- 
tion [49]. However, we suggest that the results generated 
from several independent reanalysis data sets should be 
compared if definite climatological conclusions are to be 
drawn in future studies using coupled climate subnetwork 
analysis. 

Additional information may be extracted from the avail- 
able coupled climate subnetworks by investigating the spa- 
tially fully resolved fields of local cross-network measures 
or relying on further measures such as the local cross- 
clustering coefficient (which was not shown here for brevity) . 
In a next step, network null models with a randomised 
interaction topology could be developed and used to as- 
sess the statistical significance of observed local and global 
cross-network measures. The simplest meaningful network 
null model of this type that was already discussed in the 
context of (local) cross-clustering and cross-transitivity 
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(Sec. 2.2) can be constructed by fully randomising the in- 



teraction structure between two subnetworks while keep- 
ing the number of cross-edges fixed, e.g., by first delet- 
ing all cross-edges and subsequently redistributing them 
randomly between the two subnetworks. Furthermore, the 
network construction methodology may be fine-tuned, e.g., 
by using measures for detecting nonlinear or even direc- 
tional interrelationships between time series or alterna- 
tively by including edges in the network based on the sta- 
tistical significance of their associated correlation strengths. 

Summarising the results of this first application of our 
framework for interacting network analysis, the particu- 
lar advantage of this approach is that substantial con- 
clusions can be drawn by analysing the dynamical cor- 
relation structure of the three-dimensional gcopotcntial 
height field alone, without considering fields of temper- 
ature, moisture content or other relevant climatological 
variables. Subject to future work is the application of the 
interacting networks approach to fields of distinct clima- 
tological observables (e.g., surface air temperature and sea 
surface salinity) to further investigate the coupled dynam- 
ical behaviour of different components of the climate sys- 
tem. 



4 Conclusions 



In summary, we have developed a novel graph-theoretical 
framework for investigating in detail the interaction topol- 
ogy between pairs of subnetworks embedded within a net- 
work of networks. Applying this framework to analyse 
the correlation structure of a four-dimensional (spatio- 
temporal) data set of the climatological variable geopoten- 
tial height yielded a consistent picture of the large scale 
circulation of the Earth's atmosphere. Particularly, the 
new measure cross-betweenness centrality was shown to 
have the potential to reveal previously unknown features 
of and to help address open questions on the atmosphere's 
general circulation [28] , particularly when considering its 
response to climate change [26]. Our results suggest that 
the coupled climate subnetwork approach presented in this 
work opens promising perspectives for the integrated anal- 
ysis of several fields of climatological observables or, more 
generally, spatially embedded fields of arbitrary time se- 
ries in the context of Earth system analysis. Particularly 
it will serve researchers as a tool complementary to estab- 
lished linear methods for the joint analysis of several cli- 
mate data sets like canonical correlation analysis or singu- 
lar value decomposition of the covariance matrix between 
two fields [42]. Furthermore, we expect the proposed gen- 
eral graph-theoretical framework to meet the increasing 
need for investigating and understanding the interaction 
of complex systems from domains as diverse as social sci- 
ence, technology and engineering or the life sciences, as 
well as to stimulate further research in this direction, e.g., 
into supplementing the analysis with sophisticated null 
models designed for assessing selected aspects of interact- 
ing network structure. 
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